global option k_input beta mu_z sig_z c pro z0 rho lambda alp m_bar eta

beta = 1/0.0162; % semi-demand elasticity
mu_z = 0.077;
sig_z = 0.087;
c = 0.12 ; % entry cost
pro = 8; % firm's productivity
z0 = -3.47;
rho = 0.07; % discount rate
lambda = 0.044; % exit rate
alp = 0.3; % capital share in the production function
m_bar = 1 ; % entry rate
eta = 100;
option = 'proportional'; % proportional fixed

%% baseline
k_input = [0;0;0]./100;
run 'run_aax_bunching_20210527.m'
save pre_regulation.mat  
k_input = [0.406; 0.106;0]./100;
run 'run_aax_bunching_20210527.m'
save regulation.mat 
k_input = [0.219; 0.003;0]./100;
run 'run_aax_bunching_20210527.m'
save relief.mat 

%% robustness
k_input = [0.406; 0.106;0]./100;
rho = 0.07/1.1; % discount rate
run 'run_aax_bunching_20210527.m'
save cost_of_capital.mat  
rho = 0.07; % discount rate

c = 0.12*1.1 ; % entry cost
run 'run_aax_bunching_20210527.m'
save entry_cost.mat  
c = 0.12 ; % entry cost


lambda = 0.044*1.1; % exit rate
run 'run_aax_bunching_20210527.m'
save exit_rate.mat  
lambda = 0.044; % exit rate

clearvars; 
clc; close all;

%%
dt_pre = load('pre_regulation.mat');
dt_reg = load('regulation.mat');
dt_rel = load('relief.mat');


var = {'mass','mb_all','K','R','welf_br','pi_l','pi_m','pi_h','asset_share_l','asset_share_m','asset_share_h',...
    'num_share_l','num_share_m','num_share_h'};
name = {'Mass of banks','Market-to-book','Lending quantity','Lending rate', 'Output','Small banks','Medium banks','Big banks',...
    'Small banks','Medium banks','Big banks','Small banks','Medium banks','Big banks'};
value=[];
for k=1:numel(var)
        value(k,1)=dt_pre.(var{k});
       % value(k,2)=dt_reg.(var{k});
        value(k,2)=(dt_reg.(var{k})/dt_pre.(var{k})-1)*100 ;
       % value(k,4)=dt_rel.(var{k});
        value(k,3)=(dt_rel.(var{k})/dt_pre.(var{k})-1)*100 ;
%         value(1,1)=1; % bank mass
end


   

%%
% FID = fopen('D:\Dropbox\YiyuanZhang\Counterfactuals_20200412.tex','w');
FID = fopen('D:\Dropbox\Regulation Cost\Watch what they do not what they say\FinalTables_firstdraft\Counterfactuals_20200412.tex','w');
fprintf(FID, '\\begin{centering}\\begin{tabular*}{\\linewidth}{@{\\extracolsep{\\fill}}lccc@{}}\\toprule \n');
fprintf(FID, 'Counterfactual&');
fprintf(FID, 'Baseline   & \\multicolumn{1}{c}{Dodd-Frank}    & \\multicolumn{1}{c}{Regulatory relief}     \\\\  \\midrule \n');
fprintf(FID, '  \\\\ \n');


for k = 1:5
fprintf(FID, name{k});
fprintf(FID, '&%5.3f     & %5.3f \\%%   &%5.3f \\%%   \\\\ \n',value(k,:));
end


fprintf(FID, '  \\\\ \n');
fprintf(FID, ' & \\multicolumn{3}{c}{\\textit{Annual profits} } \\\\  \\midrule \n');

for k = 6:8
fprintf(FID, name{k});
fprintf(FID, '&%5.3f     & %5.3f \\%%   &%5.3f \\%%   \\\\ \n',value(k,:));
end

fprintf(FID, '  \\\\ \n');
fprintf(FID, ' & \\multicolumn{3}{c}{\\textit{Share of assets} } \\\\  \\midrule \n');

for k = 9:11
fprintf(FID, name{k});
fprintf(FID, '&%5.3f     & %5.3f \\%%   &%5.3f \\%%   \\\\ \n',value(k,:));
end

fprintf(FID, '  \\\\ \n');
fprintf(FID, ' & \\multicolumn{3}{c}{\\textit{Share of banks} } \\\\  \\midrule \n');

for k = 12:14
fprintf(FID, name{k});
fprintf(FID, '&%5.3f     & %5.3f \\%%   &%5.3f \\%%   \\\\ \n',value(k,:));
end

fprintf(FID, '\\bottomrule \\end{tabular*} \\par\\end{centering}\n');
fclose(FID); % close tex file


%%
% Estimation result
FID = fopen('D:\Dropbox\Regulation Cost\Watch what they do not what they say\FinalTables_firstdraft\parameters.tex','w');
fprintf(FID, '\\begin{centering}\\begin{tabular*}{\\linewidth}{@{\\extracolsep{\\fill}}lllll@{}}\\toprule \n');
fprintf(FID, 'Parameter &Value & Definition   \\\\ \\midrule \n');
fprintf(FID, '$\\mu_z$ &%5.3f  & Productivity growth  \\\\  \n',dt_reg.mu_z*100);
fprintf(FID, '$\\sigma_z$ &%5.3f & Productivity diffusion   \\\\  \n',dt_reg.sig_z*100);
fprintf(FID, '$\\theta$ &%5.3f & Elasticity of funding supply  \\\\  \n',dt_reg.beta);
fprintf(FID, '$\\rho$ &%5.3f & Discount rate   \\\\  \n',dt_reg.rho*100);
fprintf(FID, '$\\lambda$ &%5.3f & Exit rate  \\\\  \n',dt_reg.lambda*100);
fprintf(FID, '$z_n$ &%5.3f & Productivity of new entrants  \\\\  \n',dt_reg.z0);
fprintf(FID, '$c_e$ &%5.3f & Entry costs  \\\\  \n',dt_reg.c);
fprintf(FID, '$A$ &%5.3f & Total factor productivity  \\\\  \n',dt_reg.pro);
%fprintf(FID, '$\\overline{m}$ &%5.3f & Mass of entry in the long-run  \\\\  \n',dt_reg.m_bar);
fprintf(FID, '$\\alpha$ &%5.3f & Curvature of the production function  \\\\  \n',dt_reg.alp);
fprintf(FID, '$\\eta$ &%5.3f & Elasticit of entry  \\\\  \n',dt_reg.eta);
% fprintf(FID, '$A$ &%5.3f & Productivity \\\\  \n',dt_reg.pro);
fprintf(FID, '\\bottomrule \\end{tabular*} \\par\\end{centering}\n');
fclose(FID); % close tex file



 %%
Q_pre = exp(dt_pre.qstar);
Q_reg = exp(dt_reg.qstar);
amin = log(3);
amax = log(250);
N = 100;
x1 = linspace(amin,amax,N)';
dx1 = x1(2) - x1(1);

err =  normrnd(0,0.04,[dt_reg.I,1]);
[bch1,dist_model_reg] = gp_dist((dt_reg.qstar + err),dt_reg.f,x1);
[~,dist_model_pre] = gp_dist( (dt_pre.qstar + err),dt_reg.f,x1);

%data_orig           = readtable('bhc_and_commercial_hh_min10_branches_for_matlab_nobranchconstraints.csv'); 
% save data_orig data_orig
load data_orig
data_orig.z         = data_orig.assets./data_orig.brnumber_hh; % Create additional variables

% Starting with correct values: -------------------------------------------
cond        = data_orig.assets >= exp(amin) & data_orig.assets <= exp(amax) & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;
fpost = ones(length(apost),1);
fpost = ones(length(apost),1);
fpre = ones(length(apre),1);
[~,dist_data_reg] = gp_dist(log(apost),fpost,x1);
[~,dist_data_pre] = gp_dist(log(apre),fpre,x1);



figure
hold on
plot(exp(bch1),(dist_data_reg./sum(dist_data_reg)),'+','LineWidth',1)
plot(exp(bch1),(dist_model_reg./sum(dist_model_reg)),'LineWidth',1)

xlabel('Asset Size')
xlim([exp(amin) exp(amax)])
ylim([min((dist_data_reg./sum(dist_data_reg))) max(dist_data_reg./sum(dist_data_reg))])
 set(gca, 'yScale', 'log')
  set(gca, 'xScale', 'log')
y = ylim; % current y-axis limits
plot(exp(dt_reg.zl(1) + dt_reg.beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
plot(exp(dt_reg.zh(1) + dt_reg.beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
plot(exp(dt_reg.zl(2) + dt_reg.beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
plot(exp(dt_reg.zh(2) + dt_reg.beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);

legend('Data','Model')
hold off

Figures = 'D:\Dropbox\Regulation Cost\Watch what they do not what they say\FinalFigures_firstdraft\';

fig = gcf;
fig.PaperPositionMode = 'manual'; % fig.PaperPositionMode = 'auto' 'auto'
fig.PaperPosition     = [1.3333    3.3125    5.8333    4.3750];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
print(fig,strcat(Figures,'/model_vs_data_3b_250b_post'),'-dpdf','-fillpage')

[b,bint,~,~,stats] =regress(log(dist_data_pre),[ones(size(bch1)),bch1])



figure
hold on
plot(exp(bch1),(dist_data_pre./sum(dist_data_pre)),'+','LineWidth',1)
plot(exp(bch1),(dist_model_pre./sum(dist_model_pre)),'LineWidth',1)

xlabel('Asset Size')
xlim([exp(amin) exp(amax)])
ylim([min((dist_data_reg./sum(dist_data_reg))) max(dist_data_reg./sum(dist_data_reg))])
 set(gca, 'yScale', 'log')
  set(gca, 'xScale', 'log')
y = ylim; % current y-axis limits
% plot(exp(dt_reg.zl(1) + beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
% plot(exp(dt_reg.zh(1) + beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
% plot(exp(dt_reg.zl(2) + beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
% plot(exp(dt_reg.zh(2) + beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
% 
legend('Data','Model')
hold off

Figures = 'D:\Dropbox\Regulation Cost\Watch what they do not what they say\FinalFigures_firstdraft\';

fig = gcf;
fig.PaperPositionMode = 'manual'; % fig.PaperPositionMode = 'auto' 'auto'
fig.PaperPosition     = [1.3333    3.3125    5.8333    4.3750];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
print(fig,strcat(Figures,'/model_vs_data_3b_250b_pre'),'-dpdf','-fillpage')



[b,bint,~,~,stats] =regress(log(dist_data_pre),[ones(size(bch1)),bch1])



%%  cdf
for i=1:6
Q_pre = exp(dt_pre.qstar);
Q_reg = exp(dt_reg.qstar);
if i<=2
    amin = i*10-1;
    amax =  i*10+1;
else
    amin = i*10-2;
    amax =  i*10+2;
end

N = 1000;
x1 = linspace(amin,amax,N)';
dx1 = x1(2) - x1(1);

err =  normrnd(0,0.04,[dt_reg.I,1]);
[bch1,dist_model_reg] = gp_dist(exp(dt_reg.qstar + err),dt_reg.f,x1);
[~,dist_model_pre] = gp_dist( exp(dt_pre.qstar + err),dt_reg.f,x1);

%data_orig           = readtable('bhc_and_commercial_hh_min10_branches_for_matlab_nobranchconstraints.csv'); 
data_orig.z         = data_orig.assets./data_orig.brnumber_hh; % Create additional variables

% Starting with correct values: -------------------------------------------
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;
fpre = ones(length(apre),1);
fpost = ones(length(apost),1);
fpre = ones(length(apre),1);
[~,dist_data_reg] = gp_dist(apost,fpost,x1);
[~,dist_data_pre] = gp_dist(apre,fpre,x1);


figure
hold on
plot(bch1,cumsum(dist_data_pre./sum(dist_data_pre)),'--','LineWidth',2)
plot(bch1,cumsum(dist_data_reg./sum(dist_data_reg)),'-','LineWidth',2)
plot(bch1,cumsum(ones(size(bch1)))./sum(ones(size(bch1))),':','LineWidth',2)
xlabel('Assets')
xlim([amin amax])
%  set(gca, 'yScale', 'log')
  set(gca, 'xScale', 'log')
  ylim([0,1])
y = ylim; % current y-axis limits
plot(exp(dt_reg.zl(1) + dt_reg.beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
plot(exp(dt_reg.zl(2) + dt_reg.beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
% legend('Post Dodd-Frank','Pre Dodd-Frank','Power law')
legend('Pre Dodd-Frank','Post Dodd-Frank','Location', 'southeast')
hold off

fig = gcf;
fig.PaperPositionMode = 'manual'; % fig.PaperPositionMode = 'auto' 'auto'
fig.PaperPosition     = [1.3333    3.3125    5.8333    4.3750];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
print(fig,strcat(Figures,'/cdf_data',num2str(i)),'-dpdf','-fillpage')
end



%%  cdf post relief
for i=1:6
Q_pre = exp(dt_pre.qstar);
Q_reg = exp(dt_reg.qstar);
if i<=2
    amin = i*10-1;
    amax =  i*10+1;
else
    amin = i*10-2;
    amax =  i*10+2;
end

N = 1000;
x1 = linspace(amin,amax,N)';
dx1 = x1(2) - x1(1);

err =  normrnd(0,0.04,[dt_reg.I,1]);
[bch1,dist_model_reg] = gp_dist(exp(dt_reg.qstar + err),dt_reg.f,x1);
[~,dist_model_pre] = gp_dist( exp(dt_pre.qstar + err),dt_reg.f,x1);

%data_orig           = readtable('bhc_and_commercial_hh_min10_branches_for_matlab_nobranchconstraints.csv'); 
data_orig.z         = data_orig.assets./data_orig.brnumber_hh; % Create additional variables

% Starting with correct values: -------------------------------------------
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 2 | data_orig.post ==3); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.post == 2,:).assets;
apost = data(data.post == 3,:).assets;
fpre = ones(length(apre),1);
fpost = ones(length(apost),1);
fpre = ones(length(apre),1);
[~,dist_data_reg] = gp_dist(apost,fpost,x1);
[~,dist_data_pre] = gp_dist(apre,fpre,x1);


figure
hold on
plot(bch1,cumsum(dist_data_pre./sum(dist_data_pre)),'--','LineWidth',2)
plot(bch1,cumsum(dist_data_reg./sum(dist_data_reg)),'-','LineWidth',2)
plot(bch1,cumsum(ones(size(bch1)))./sum(ones(size(bch1))),':','LineWidth',2)
xlabel('Assets')
xlim([amin amax])
%  set(gca, 'yScale', 'log')
  set(gca, 'xScale', 'log')
  ylim([0,1])
y = ylim; % current y-axis limits
plot(exp(dt_reg.zl(1) + dt_reg.beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
plot(exp(dt_reg.zl(2) + dt_reg.beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
% legend('Post Dodd-Frank','Pre Dodd-Frank','Power law')
legend('Pre relief','Post relief','Location', 'southeast')
hold off

fig = gcf;
fig.PaperPositionMode = 'manual'; % fig.PaperPositionMode = 'auto' 'auto'
fig.PaperPosition     = [1.3333    3.3125    5.8333    4.3750];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
print(fig,strcat(Figures,'/cdf_data_relief_',num2str(i)),'-dpdf','-fillpage')
end


%%  pdf
for i=1
Q_pre = exp(dt_pre.qstar);
Q_reg = exp(dt_reg.qstar);
if i<=2
    amin = i*10-5;
    amax =  i*10+5;
else
    amin = i*10-2;
    amax =  i*10+2;
end

N = 30;
x1 = linspace(amin,amax,N)';
dx1 = x1(2) - x1(1);

err =  normrnd(0,0.04,[dt_reg.I,1]);
[bch1,dist_model_reg] = gp_dist(exp(dt_reg.qstar + err),dt_reg.f,x1);
[~,dist_model_pre] = gp_dist( exp(dt_pre.qstar + err),dt_reg.f,x1);

%data_orig           = readtable('bhc_and_commercial_hh_min10_branches_for_matlab_nobranchconstraints.csv'); 
data_orig.z         = data_orig.assets./data_orig.brnumber_hh; % Create additional variables

% Starting with correct values: -------------------------------------------
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;
fpre = ones(length(apre),1);
fpost = ones(length(apost),1);
fpre = ones(length(apre),1);
[~,dist_data_reg] = gp_dist(apost,fpost,x1);
[~,dist_data_pre] = gp_dist(apre,fpre,x1);


figure
hold on
plot(bch1,(dist_data_pre./sum(dist_data_pre)),'+','LineWidth',2)
plot(bch1,(dist_data_reg./sum(dist_data_reg)),'-','LineWidth',2)
% plot(bch1,(ones(size(bch1)))./sum(ones(size(bch1))),':','LineWidth',2)
xlabel('Assets')
xlim([amin amax])
%  set(gca, 'yScale', 'log')
  set(gca, 'xScale', 'log')
%   ylim([0,1])
y = ylim; % current y-axis limits
plot(exp(dt_reg.zl(1) + dt_reg.beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
plot(exp(dt_reg.zl(2) + dt_reg.beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
% legend('Post Dodd-Frank','Pre Dodd-Frank','Power law')
legend('Pre Dodd-Frank','Post Dodd-Frank','Location', 'southeast')
hold off

fig = gcf;
fig.PaperPositionMode = 'manual'; % fig.PaperPositionMode = 'auto' 'auto'
fig.PaperPosition     = [1.3333    3.3125    5.8333    4.3750];
fig_pos = fig.PaperPosition;
fig.PaperSize = [fig_pos(3) fig_pos(4)];
% print(fig,strcat(Figures,'/cdf_data',num2str(i)),'-dpdf','-fillpage')
end

